# nohup /usr/local/R/2.15.3/bin/R CMD BATCH HR_treaty04.R &


# set time start variable
systime1 <- Sys.time()
print(Sys.time() - systime1)


# set time start variable
time1 <- Sys.time()
print(Sys.time() - time1)

# load libraries
#library(rjags); library(coda); library(foreign)
.libPaths( c( .libPaths(), "/usr/lib64/R/library") )
library(rjags); library(coda); library(foreign)


#----------------------------------------------------------------------------------------------------#

data <- read.csv("hr_treaties.csv")
data <- subset(data, year>=1948)

# subset out international humanitarian law instruments and keep only international human rights law instruments
data <- subset(data, select=c(-gc49_rn, -cslwar_rn, -gcop177_rn, -gcap277_rn, -catop_rn))

# subset out regional human rights law instruments
data <- subset(data, select=c(-echr_rn, -iachr_rn, -achpr_rn, -iacat_rn, -ecat_rn))



year <- NA
year[1] <- 1
country <- NA
panel <- NA
panel.count <-1
i <- 2
country[1] <- 1
j <- 1
while(i <= nrow(data)){
if(data$ccode[i]!=data$ccode[i-1]){
panel[j] <- panel.count
panel.count <- 0
j <- j+1
}
country[i] <- j
#i <- i + 1
panel.count <- panel.count + 1
year[i] <- panel.count
i <- i + 1
}
j
panel[j] <- panel.count

# full model
#y.pre <- as.matrix(data[,4:36])
#y.pre <- as.matrix(data[,4:31])
y.pre <- as.matrix(data[,4:ncol(data)])
y <- y.pre

head(y)

n <- nrow(y.pre)
n

k <- ncol(y.pre)
k


# random initial values
MakeInits <- function(){
 # nIRT.Binary <- length(4:36)
 # nIRT.Binary <- length(4:31)
  nIRT.Binary <- length(4:ncol(data))
  MU <- matrix(rnorm(length(panel)*max(panel), mean = 0, sd = 1), nrow=length(panel),  ncol=max(panel))

  BETA <- runif(nIRT.Binary)
  ALPHA <- runif(nIRT.Binary)

  SIGMA <- runif(1)

  out <- list(mu=MU, alpha=ALPHA, beta=BETA, sigma=SIGMA)
  return(out)
}
inits1 <- MakeInits()
inits2 <- MakeInits()
inits3 <- MakeInits()

inits.function <- function(chain){
  return(switch(chain,
         "1"=inits1,
         "2"=inits2,
         "3"=inits3
         ))
}

# jags.model arguments
ADAPT <- 1000
BURNIN <- 5000
DRAWS <- 15000
THIN <- 15
CHAINS <- 3
NAME <- "HR_treaty04.bug"

# print time taken
print(Sys.time() - time1)
  
#rjags code version
m <- jags.model(file=NAME, data=list("y"=y, "year"=year, "country"=country, "n.country"=length(panel), "n.year"=max(panel), "n"=nrow(y), "k"=k), inits=inits.function, n.chains=CHAINS, n.adapt=ADAPT)

print(Sys.time() - time1)



time1 <- Sys.time()
print(Sys.time() - time1)


update(m, BURNIN)
print(Sys.time() - time1)


# take posterior draws
time1 <- Sys.time()
print(Sys.time() - time1)

M <- coda.samples(m, variable.names=c("x", "alpha", "beta", "kappa", "sigma"), n.iter=DRAWS, progress.bar="text", thin=10)


print(Sys.time() - time1)


save.image(file="image.Rdata")


print(Sys.time() - systime1)



mat1 <- as.matrix(as.mcmc(M[[1]]))
mat2 <- as.matrix(as.mcmc(M[[2]]))
mat3 <- as.matrix(as.mcmc(M[[3]]))
posterior.estimates <- rbind(mat1, mat2, mat3)
write.csv(as.data.frame(posterior.estimates), "EstimateLatentHRtreaty03.csv", row.names=F)

parameter.mean <- apply(posterior.estimates, 2, mean)
parameter.sd <- apply(posterior.estimates, 2, sd)

save.image(file="image.Rdata")


# make country year file
data$latentmean <- parameter.mean[(length(parameter.mean)-nrow(y)+1):length(parameter.mean)]
data$latentsd <- parameter.sd[(length(parameter.sd)-nrow(y)+1):length(parameter.sd)]
write.csv(data, "CYEstimateLatentHRtreaty04.csv", row.names=F)


save.image(file="image.Rdata")


